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ABSTRACT 


A comparison  is  made  between  several  different  methods  for 
calculating  energy  transport  within  a wave  field.  Two  Langevin 
techniques  are  developed.  The  first  is  based  on  the  fluctuation- 


dissipation  theorem  and  provides  relaxation  rates  and  a 
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transport  equation.  The  second  method  is  an  application  of 
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the  Krylov-Bogoliubov-Mitropolsky  perturbation  theory  and  pro- 
vides a Langevin  rate  constant  Vp  at  lowest  order.  The  two 
formulations  are  shown  to  be  closely  related  to  the  radiative 
transfer  (Boltzmann)  equation  whose  rate  is  the  difference 
between  ftp  and  Vp.  Specific  application  of  the  Langevin 
methods  is  to  internal  waves  in  the  ocean.  Computations  show 
that  the  GM-76  spectrum  is  approximately  an  equilibrium 
spectrum  except  for  frequencies  near  the  inertial  frequency  and 
at  the  lowest  vertical  modenumbers.  The  sensitivity  of  Vp 
and  Vpto  spectral  form  is  also  discussed.  Simple  analytic 
expressions  for  the  rates  are  derived  for  the  induced  diffusion, 
elastic  scattering,  and  parametric  subharmonic  instability 
mechanisms.  Only  the  first  of  these  mechanisms  is  ever  of  much 
numerical  significance.  Finally,  net  energy  flow  in  thg.  non- 
equilibrium portion  of  the  GM-76  spectrum  is  discussed. 
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1.  INTRODUCTION 


Nonlinear  energy  transfer  mechanisms  within  the  oceanic  internal 
wave  field  have  been  studied  by  Olbers  [1976]  and  MoComas  and  Bretherton 
[1977].  These  authors  used  a radiative  transfer  equation  [ Hasselmann , 
1966,  1967]  for  their  computations.  The  radiative  transfer  (or  Boltzmann) 
equation  governs  the  evolution  of  wave  action  spectra  ensemble  averaged 
over  many  realizations  of  the  wave  field.  Derivation  of  this  equation 
from  the  fluid  equations  requires  several  approximations: 


(a)  Nonlinearities  are  assumed  "weak,"  and  only  lowest  order 
(quadratic  for  the  case  of  internal  waves)  nonlinear  terms  are  retained 
in  the  equations  of  motion. 


(b)  Two-time  perturbation  methods  are  used;  the  "fast  time" 
corresponds  to  linear  wave  periods  and  the  "slow  time"  to  nonlinear 
interaction  timescales. 


(c)  Spatial  homogeneity  in  any  horizontal  plane  is  assumed.  This 
allows  simplification  of  second  moments  of  wave  amplitudes. 


(d)  A closure  approximation  is  made  by  the  discard  of  fourth  and 
higher  order  cumulants. 


In  the  present  study  we  compare  several  different  techniques  for 


calculating  energy  transport  within  the  internal  wave  field.  The  dynamical 
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equations  [Meiss,  Pomphrey , and  Watson ; 1979,  henceforth  referred  to  as  I] 
are  derived  using  approximation  (a)  and  describe  nonlinear  interaction 
of  the  linear  normal  modes  of  the  wave  field.  Here  we  use  the  Garrett- 
Munk  exponential  profile  for  the  Vaisala  frequency  and  the  WKB  approxi- 
mation to  calculate  the  vertical  modefunctions  {Others  [1976]  and 
McComas  and  Bretherton  [1977]  have  used  a constant  Vaisala  profile). 

Two  Langevin  equation  techniques  are  developed  in  this  paper  to  study  the 
dynamics  of  the  system  derived  in  I (Langevin  methods  have  been  used  by 
Holloway  and  Hendershott  [1977]  to  discuss  Rossby  waves). 

The  first  Langevin  method  is  based  on  the  fluctuation-dissipation 
theorem  [Lax,  I960,  1966].  This  method  provides  relaxation  rates  and 
a transport  equation,  and  depends  rather  little  on  dynamics.  It  assumes 
linear  relaxation  to  a known  equilibrium  state  and  also  requires  approxi- 
mations (b),  (c)  and  (d).  For  the  "known  equilbrium  state,"  we  choose 
a set  of  Garrett-Munk  related  spectra,  with  the  "GM-76"  spectrum 
[ Cairns  and  Williams,  1976]  as  our  reference  standard. 

To  obtain  the  second  form  of  the  Langevin  equation  we  use  approxi- 
mation (b)  in  the  form  of  the  Krylov-Bogoliubov-Mitropolsky  perturbation 
method.  [ Bogoliubov  and  Mitropolsky,  196l].  Approximation  (d)  is  not 
required  to  calculate  the  Langevin  rate  constant  with  this  procedure. 

We  shall  show  in  Section  3 that  the  two  Langevin  formulations  are 
closely  related  to  the  radiative  transfer  equation.  For  fluctuations  near 
a true  equilibrium  state  the  two  formulations  are  in  agreement.  In  this 
case  the  transport  equation  obtained  from  the  Langevin  equation  is 
identical  in  form  with  that  used  by  Others  [1976]  and  by  McComas  and 
Bretherton  [ 1977  ] . 


I 
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The  Langevin  method  leads  to  a decay  rate  v(k)  for  the  autocorrelation 
of  the  amplitude  for  a linear  internal  wave  labelled  as  "k."  When  calcu- 
lated using  the  fluctuation-dissipation  theorem,  we  call  this  Vp(k). 

The  value  obtained  from  perturbation  theory  is  written  as  Vp(k).  The 
radiative  transfer  equation  used  by  Olbevs  [1976]  and  by  MoComas  and 
Bretherton  [1977]  for  the  action  density  <J  > is  of  the  form 

" it 

« <Jk>  = 2 VB(k)  "V  ’ 


where 

Vg(k)  = Vp(k)  - Vp(k) 

It  is  convenient  to  think  of  as  representing  the  rate  of  energy  input 
to  mode  k from  the  "noise  field"  of  the  wave  system  and  to  think  of 
Vp  as  describing  a rate  of  energy  loss  from  mode  k. 

MaComas  [1978]  reported  numerical  experiments  in  which  he  intro- 
duced small  distortions  in  the  GM-76  spectrum  and  computed  the  relaxation 
to  "equilibrium."  We  give  an  analytic  description  of  this  in  Section  3 
and  show  that  the  Langevin  autocorrelation  decay  rate  vp(k)  also  determines 
the  rate  of  return  to  equilibrium. 

Numerical  calculations  of  Vp  and  Vp  are  presented  in  Section  5.  We 
shall  see  that  for  the  GM-76  spectrum,  except  for  frequencies  near  the 
inertial  frequency  and  at  the  lowest  vertical  mode  numbers,  |vg|  is 
from  one  to  three  orders  of  magnitude  smaller  than  Vp  (or  Vp) . We 
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thus  find  precise  numerical  confirmation  of  the  conclusion  of 
McComas  and  Bretherton  [1977]  (who  calculate  only  vp)  that  GM-76 
is  approximately  an  equilibrium  spectrum. 

The  sensitivity  of  and  to  spectral  form  is  also  discussed 
in  Section  5 for  a class  of  GM-76  type  spectra.  Not  surprisingly , 
in  the  eauilibrium  region  where  v_  is  a small  difference  between  two 

D 

much  larger  quantities,  the  values  of  vp  are  rather  sensitive  to  spectral 
shape . 

McComas  and  Bretherton  [1977]  describe  three  limiting  mechanisms 
for  energy  transfer.  They  call  these  induced  diffusion,  elastic  scattering 
and  parametric  subharmonic  instability.  In  Section  6 we  present  simple 
analytic  expressions  for  and  vp  for  each  of  these  mechanisms.  We 
show  from  these  expressions  that  GM-76  represents  an  equilibrium  spectrum, 
except  at  the  lowest  frequencies  and  lowest  vertical  mode  numbers , with 
respect  to  both  induced  diffusion  and  elastic  scattering.  (For  the  case 
of  induced  diffusion  this  was  noted  by  McComas  and  Bretherton  [1977]  and 
was  implied  more  generally  from  their  numerical  calculations.) 

We  show  that  the  elastic  scattering  and  parametric  subharmonic 

instability  mechanisms  are  never  of  much  numerical  significance  for  the 

GM-76  spectrum.  For  frequencies  greater  than  three  times  the  inertial 

frequency  the  analytic  expressions  obtained  from  induced  diffusion  for 

both  and  v_  agree  well  with  the  numerically  calculated  values  for 
F P 

these  quantities. 

Net  energy  flow  in  the  non-equilibrium  portion  of  the  GM-76  spectrum 


is  discussed  in  Section  7-  Our  conclusions  are  generally  consistent  with 
those  of  Others  [1976]  and  of  McComas  and  Bretherton  [1977]  that  energy 
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is  transferred,  from  the  low  vertical  mode  number  regime  into  that  of  high 
vertical  mode  numbers  and  near  inertial  frequencies.  Since  net  energy- 
flow  is  determined  by  v , the  details  depend  relatively  sensitively  on 

■D 

the  assumed  spectrum. 

The  Langevin  rate  constants  discussed  here  govern  relaxation  processes 
in  the  wave  field  and  also  immediately  yield  values  for  the  more  frequently 
calculated  Boltzmann  constant  The  simple  analytic  approximations  for 

these  rate  constants  given  in  this  paper  should  be  adequate  for  applications 
to  internal  wave  transport  processes.  j- 


, 


I 

I 
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2.  THE  DYNAMICAL  MODEL 

A general  description  of  internal  wave  phenomena  and  theory  is 
given  by  Phillips  [1977].  The  detailed  description  of  the  specific 
dynamical  model  used  here  is  given  in  I.  A rectangular  coordinate 
system  is  chosen  with  the  x-y  plane  tangent  (locally)  to  the  ocean 
surface.  The  bottom  is  assumed  horizontal  at  z = -K.  The  Garrett-Munk 
exponential  V&is&ia  profile  is  used  in  this  paper: 

N(z)  = Nq  exp(z/B)  . (1) 

Dimensional  quantities  are  N = 5.2  (10)  rad/sec,  B = 1.2  km  and  the 

o 

surface  fluid  density  p = p(o).  The  coriolis  frequency  is  assumed 

o 

vertical  with  magnitude 

(2) 

f = 7.3  (10)"5  rad/sec  = .OlU  II 

0 

corresponding  to  30°  latitude.  Vectors,  x = (x,y),  are  two-dimensional 
in  the  horizontal  plane. 

The  vertical  displacement  of  a Lagrangian  fluid  element  was  written 
in  I as 


C3(x,z,t)  = Re  [Y(x,z ,t ) ] , 


(t) 


Y(x,z,t)=iN  /B  7 y W,  (z)  e 

o Yi  r k ka 

a=l  k 


ik*x 


(3) 


Equation  (3)  is  a Fourier  expansion  in  an  ocean  of  rectangular  area 


o 
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in  terms  of  horizontal  wavenumbers  k.  The  linear  vertical  modefunctions , 
Wjw/z),  satisfy  eigenval^e  equation 


f ! p f W ] + pk  0 -r  W = 0, 
dz  dz  ka  2,,  N „2  ka 

aj  (kj-f 
a 


N2-w2(k) 


to  (k)>0, 
a 


W,  (0)  = V ( -H)  = 0. 
ka  ka 


Here  u>a(k)  is  the  angular  frequency  of  a linear  internal  wave  with 
integer  vertical  modenumber  a. 

The  amplitudes  a^  are  dimensionless  wave  slope  variables  and  satisfy 
[I,  eq.  2.23] 


a,  + iw  a,  = y [6,  G^a*a  +5  G^  a a ] 
K alt  „L  k+H-m  m l m k-H-m  to  l m 

Q m 


Here  we  have  used  abbreviated  labels,  writing  k,!,m  for  (k,a), 
and  (m,y)  respectively.  Explicit  expressions  for  the  "G"  coupling 
coefficients  in  (5)  were  given  in  I.  They  contain  integrals  over  the 
product  of  three  modefunctions  W (see  appendix  A)  as  well  as  factors 
determined  by  the  geometry  of  the  interacting  triad. 

Wave  action  per  unit  area  is  expressed  in  terms  of  the  slope 
variables  as  [I,  eq.  2.21] 


, (o  S '^a1 

J,  = p N B3  ” 0 — — - . 

00  w^-f2  2(kB)2 

a 
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The  Boussinesq  approximation  has  not  been  used  to  derive  equations 
(3),  (1+)  and  (5).  This  approximation  entails  the  neglect  of  the  term 

d dWka 

in  (4).  In  Appendix  A it  is  shown  to  have  a negligible 
effect  on  the  coupling  coefficients  in  (5).  We  shall  therefore  make 
this  simplifying  approximation  for  our  subsequent  discussion. 

Following  Garrett  and  Munk  [1975,  1979]  we  use  the  WKB  approximation 
to  solve  (4).  This  appears  reasonably  valid  for  a>2  (results  for  a=l 
are  at  best  qualitative  but  will  be  included  for  completeness).  We  also 
use  their  dispersion  relation 


w2(k)  = f2  + 


(a-Jfi)ir 


valid  when  u (k)<<N  (and  assumed  valid  in  this  paper  for  u>  (k)<N  /3). 

a o a o 

The  power  spectral  density  (PSD)  of  vertical  displacement  is  written  as 
<jj(k,ct,z)  and  normalized  so  that 


2 2 

<£,>  = l fd  k Ij;(k,a,z)  , 


where  "<  >"  represents  an  ensemble  average  over  realizations  of  the  internal 

wave  field.  We  shall  require  only  the  spectrum  iji(k,a)  extrapolated  to 

2 2 2 

the  surface  for  which  <£(r>  = <s  > = (7- 3m)  . From  equations  (l),  (3), 

i o 

(6),  (7)  and  the  normalization  of  the  eigenfunctions  W^  (Appendix  A), 
we  see  that 


<l»(k,a)  = 


"'■kJ  " 


In  this  paper  we  adopt  a PSD  related  to  the  Garrett-Munk  form  [Garrett 
and  Munk , 1975,  1979],  with 

(kB  )a  P_1 

o o o P/2  + 1 T 

[(kB)2  + 1.9  (10)"  a ] [1  + (a/3)  ] 

(10) 

where  p is  the  "wavenumber  slope"  and  t is  the  "modenumber  slope."  The 
GM-76  spectrum  [Cairns  and  WiZZiams,  1976]  with  p=2,  t=2,  and  W ^ = .013 
is  chosen  as  our  "standard." 
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3.  RELAXATION  RATES 

In  this  section  we  discuss  and  relate  three  different  methods  for 
calculating  relaxation  rates  in  a random  wave  field.  We  fix  attention 
on  a single,  definite  internal  wave  mode,  say  that  labelled  (k,a).  This 
is  the  "test  wave,"  and  we  study  its  interaction  with  the  ambient  wave- 
field.  Averages  are  denoted  "<  >"  and  are  over  an  ensemble  of  states 
of  the  ambient  field.  We  suppose  that  the  test  wave  always  has  a definite 
amplitude  at  some  initial  time,  say  t=0.  As  t-*»  there  will  be  no  difference 
between  this  ensemble  and  an  ensemble  of  states  of  the  entire  internal 
wave  field  (including  the  test  wave). 

A.  Fluctuation  Dissipation  Theorem 

We  cast  the  equations  of  motion  (5)  for  the  test  wave  as  a Langevin 
equation  by  representing  the  nonlinear  terms  on  the  right  hand  side  by 
a random  force  R(t)  (for  convenience  we  temporarily  drop  the  mode  labels 
(k,a)) : 

a + iwa  = R(t)  . (ll) 

From  this  point  of  view  the  test  wave  is  driven  by  the  ambient  waves 
which  act  as  an  "equilibrium  heat  bath."  It  is  convenient  to  transform 
variables,  defining 
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b(t)  = eiut  a(t ) ; 


so  that  (ll)  becomes 


b = elu)t  R i R(t) 


Following  the  conventional  argument  [Lax,  i960,  1966]  we  suppose  that 
at  t=0,  b has  the  definite  value  b(0).  The  mean  of  b,  averaged  over  the 
ambient  ensemble,  obeys  the  equation 


fi  <h>  - • 


It  is  anticipated  that  as  t-**>,  <b>  -*•  0.  This  leads  to  a reasonable 
postulate  for  the  form  of  R: 


R(t)  = - vb(t)  + F(t)  , 


<F>  = 0 


Here  v is  the  Langevin  relaxation  rate  (assumed  real  since  any  imaginary 
part  may  be  removed  by  a transformation  of  the  form  ( 12 )) . 

From  assumption  (l^),  equation  (13)  becomes 


b + vb  = F(t) 


(15) 
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from  which  it  follows  that 


<b(t )>  = b(0)  e 


To  develop  the  fluctuation  dissipation  theorem  we  assume 


:|b(t). 


a constant,  equilibrium  value.  (Since  the  hamiltonian  from  which  equation 
(5)  was  derived  has  no  lower  bound  in  energy,  a true  equilibrium  in  the 
thermodynamic  sense  does  not  exist.  However,  large  amplitude  fluctuations 
are  sufficiently  rare  that  we  can  ignore  this  problem.  Note  the  analogy  with 
the  theory  of  low-lying  Stark  states  of  hydrogen.)  It  is,  of  course, 
of  interest  to  determine  if  the  Garrett-Munk  PSD  (10)  corresponds  to 
this  equilibrium,  and  this  is  one  of  the  tasks  for  the  computations  in 


Sec.  5. 


Considerable  simplification  results  by  assuming  6-correlation  of  the 


<F(t)F*(t' )>  = 2D  6 ( t-t * ) . 


Numerical  investigations  described  later  indicate  that  decorrelation  is 
sufficiently  rapid  that  equation  (l8)  is  a good  approximation  for  our 
application.  It  is  easy  to  show  that  D is  real. 


kJL 


Integration  of  (15)  using  assumptions  (IT)  and  (l8)  and  restoring 
mode  labels  yields  the  fluctuation-dissipation  result 


D(k,ot) 


v_(k,a)  = a 


(19a) 


where  the  subscript  "F"  indicates  that  this  relaxation  rate  is  derived 
by  the  fluctuation-dissipation  method. 

For  future  comparison  with  the  radiative  transfer  equation  we 
generalize  (19a)  to  define  a rate  coefficient 


Vp(k,a)  = 


D(k,a) 


(19b) 


where  <lajcaH>  does  not  necessarily  represent  an  "equilibrium"  value. 

To  evaluate  D we  use  the  linear  approximation  for  the  ambient  field 
amplitudes  on  the  right  hand  side  of  (5) 


a16(t)  ' a».S(0)  e 


-l“s* 


We  also  use  the  cumulant  discard  approximation  to  reduce  fourth  moments 
to  second  moments.  Evaluation  is  straightforward,  giving 


D(k,o)  = / <F(t)F#(0)>  dt 


= it  £ { I G ^ I 6,  0 6 ( to  +u>  -u)  ) 

nL  1 m 1 k+£-m  a 8 Y 


+ 2|G^  I2  6.  . 6 ( a)  -u)Q-u>  )}  <|a0J2>  <|a  |2>  . (20) 

1 Im1  k-Sl-m  a B y ’ 2.0 1 1 my 1 
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For  a large  ocean  area  T.  we  may  replace  wavenumber  sums  by  integrals 

o 

with  the  substitution 


I * -%  /A 

Z UlT 


and  use  equations  (20)  and  (9)  to  re-write  (19b): 


(k,a)  = 2ir  £ fd2ld2m  { | Gk ^ | 2 6(k+£-m)  5(gj  +w  -u  ) 


8,Y 


m 


a 3 Y 


2 2 

+ 2 1 G^m | 2 6(k-£-m)6Ua-a)g-u>Y)}  -Mj- 


£,8  )i|>(m,y ) 

iJ>(k,oi) 


It  is  clear  that  v„  is  a positive  quantity, 
r 


(21) 


i 


4 


3.  Perturbation  Method 

The  Krylov-Bogoliubov-Mitropolsky  [1961]  two-time  perturbation  method 
provides  an  alternative  means  of  obtaining  a Langevin  equation  of  the 
form  (15)  [see  e.g.  Case,  1966].  For  the  lowest  order  perturbation 
calculation,  the  cumulant  discard  assumption  is  not  required  since  the 
test  wave  is  a sure  quantity  at  t=0  under  the  averages.  Using  an 
adiabatic  "switch  on"  in  the  remote  past  of  the  right  hand  side  of  (5), 
we  obtain  a complex  "frequency  shift" 


■ re  ■ 


km  ink*  orX  ml 

+ G*  Ge  + 2G^mGk  . <|a  (0)|2> 

<o  -10  +u)  +iri  k-a+m  <o  -qj  -to  +in  k-Jl-m  ' £3^  ' 

a 6 y > _ „ a 6 Y - - - 


The  real  part  of  this  is  the  relaxation  rate 


vp(k,a)  = Re[v(k,a) ] 


= - 2ir  l fd2ld2m  2G^G“.  6(k+£-m)6(«o  +<oD-(o  ) 

„ m Ilk  . . . a 8 Y 

b»y 


+ G G 6(k-£+m)6(co  -to  +tu  ) 
ft  x,  _ _ _ a 8 Y 


Ir  jiifl  O 

+ 2G,  G * 6(k-*-m)6(w  -toD-to  ) l 
)tm  k , _ _ a 8 Y 


Here  we  have  used  the  standard  relation 


Re  [ lim  — — : — ] 

n+o(+)  x+lri 


tt6(x) 


and  also  expression  (9)  for  the  PSD.  In  Appendix  (B),  equation  (23)  is  re- 
written in  a form  which  shows  that  the  first  and  last  terms  are  positive 


v(k,a)  = - 


r 2Gk*G“ 

i lim  I - — — 6. 

n+0(+)  Z,m  V“B-Vln 


while  the  second  term  is  negative. 


| 
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C.  Radiative  Transfer  Equation 

The  radiative  transfer  or  Boltzmann  equation  has  been  used  by 

Others  [1976]  and  McComas  and  Bretherton  [1977]  to  describe  nonlinear 

internal  wave  interactions.  In  this  method,  the  equations  of  motion 

(5)  are  used  to  obtain  an  equation  for  d<|a^|  >/dt  [see,  for  example, 

Davidson,  1972]  in  terms  of  <akaj^a*>,  etc.  Equations  for  the  rate  of 

change  of  these  third  order  moments  involve  fourth  moments.  Closure 

results  from  discard  of  fourth  order  cumulants,  leaving  second  order 

moments.  Use  of  the  homogeneity  assumption  allows  expression  of  these 

in  terms  of  the  PSD.  Finally,  first  order  perturbation  theory  is  used 

to  integrate  the  equations  for  <a,  a„a*>.  The  result  of  all  this  is  the 

k 1 m 

transport  equation 


Iff 


where 


2vB(k,a) 


v 


B 


(24a) 

(24b) 


The  expression  (21)  and  (23)  for  Vp  and  \)p  are  to  be  used  here.  Symmetry 
properties  of  the  G-coefficients  may  be  used  to  rewrite  (24)  in  precisely 
the  form  given  by  Others  [1976]  and  McComas  and  Bretherton  [1977] 

(Appendix  B). 

Equation  (24)  shows  the  relationship  of  the  Boltzmann  transport 
equation  to  the  two  versions  of  the  Langevin  description.  We  shall  see 
in  the  next  section  that  (24)  may  also  be  derived  directly  from  the  Langevin 
equations. 


!J 
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In  the  application  of  the  fluctuation-dissipation  theorem,  we  were 

2 

required  to  assume  that  the  <lajcal > corresponded  to  the  equilibrium 

state  . This  was  not,  of  course,  used  to  obtain  equations  (23)  and 

(24),  which  result  from  dynamical  equations.  We  see  from  (24),  however, 
2 

that  if  <laiCctK>  = °ka»  "then  Vg(k,a)  = 0 and 


vp ( k , a ) = vF(k,a)  = vp(k,a) 


(25) 


The  two  forms  of  the  Langevin  equation  and  the  transport  equation  (24) 
are  then  consistent.  We  have  seen  that  the  quantity  vp  is  always  positive 
Evidently,  for  an  equilibrium  spectrum  vp  must  also  be  positive  although 
this  is  not  guaranteed  since  the  full  expression  (23)  for  vp  is  not 
positive  definite. 

When  the  ambient  field  is  at  (or  nearly  at)  equilibrium,  so  vp>0, 
the  fluctuation-dissipation  "noise"  D tends  to  excite  the  mode  (k,a) 
at  a rate  v-,  (or  vL).  Energy  is  lost  to  the  ambient  field  at  the  rate 

r r 

vp . To  illustrate  the  implications  of  this,  suppose  that  our  test  wave 
has  arbitrary  initial  amplitude,  but  all  other  modes  are  in  equilibrium. 

We  can  then  use  (19b)  to  integrate  (24): 


o P "2vpt 

<laka(t)I2>  = <laka(0)|2>  e + 


1-e 


*2vpt 


(26) 


The  equilibrium  value  is  evidently 

D(k,a) 
vp(k,a) 


^ - - 


n«~ 
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as  required  by  the  fluctuation-dissipation  theorem  and  (25).  Now, 
from  (l6),  we  obtain 

-v„t-iw  x t 

‘V*41  *L<0)>  ' 'iV*0’!  =•  9 ' “ ' (27> 

So  from  (26),  the  time  scale  for  <laka!  > to  reach  equilibrium  is 

determined  by  Vp,  while  the  autocorrelation  function  (27)  decay  rate 

is  determined  by  v . 

r 

The  dynamical  calculations  reported  in  I,  corresponding  to  numerical 
integration  of  equations  (5)  seem  to  be  reasonably  consistent  with  (26): 

A test  wave  of  initially  small  amplitude  was  found  to  grow  to  the  GM-76 
value  in  roughly  the  expected  time.  It  remained  at  this  level  as  long  as 
the  calculation  was  continued. 

In  the  above  discussion  we  assume  that  there  exists  an  equilibrium 
solution  to  (2U).  For  this  to  have  physical  interest  the  solution  should 
correspond,  to  some  extent,  to  observed  internal  wave  spectra.  In 
Section  5 we  shall  present  computed  values  of  Vp  and  Vp  for  a class  of 
GM  spectra  and  conclude  as  did  MoComas  and  Bretherton  [1977]  that 
GM-76  is  nearly  an  equilibrium  spectrum,  except  for  frequencies  close 
to  the  inertial  frequency  and  for  the  lowest  mode  numbers . Within 
the  "equilibrium  range"  we  will  see  that  v may  be  several  orders  of 

D 

magnitude  smaller  than  v , and  therefore  tends  to  be  quite  sensitive 
( to  details  of  the  spectrum. 


-20- 


U . F0KKER-PLANCK  EQUATION 

The  assumption  that  the  "noise"  term  in  (15)  represents  a 
Markoff  process,  fluctuating  rapidly  on  the  time  scale  v”'L,  permits  one 
to  derive  a Fokker-Planck  equation  for  the  probability  distribution  of 
the  amplitude  b [ Chandrasekhar , 19^3;  Wang  and  Uhlenbeak,  191+5].  To 
obtain  this,  we  first  write  b in  terms  of  its  real  and  imaginary  parts, 

b = x + i y 

The  probability  density  for  x and  y at  time  t is  written  as  P(x,y,t). 
The  Fokker-Planck  equation  is  obtained  from  (15): 


f - k [vF x + b tvpJ 


3x 


[f]  , 


(28) 


Here  D is  the  quantity  defined  by  (20).  An  equation  of  evolution  for 
the  wave  amplitude  intensity 


< | a | ^>  = f(x2  + y2)  P dxdy 


may  be  obtained  immediately  from  (28): 


<|akJ2>  = 2Vk’a)  " <la>«|2>! 


F . ka  1 ka1 


_d_ 

dt 


(29) 
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This  is  precisely  equivalent  to  (2U)  if  there  exists  an  equilibrium 

snectrum  o,  since  then  v_  = v . In  Appendix  3 (29)  is  rewritten  to 

r r 

include  the  effects  of  inhomogeneous  media. 

We  have  now  introduced  three  different  rate  constants,  v?,  v^,  and 

v_.  It  is  to  be  observed  from  equations  (2U),  (26),  (27)  and  (29)  that 
3 

v_  describes  the  net  rate  of  transfer  of  action  (or  energy),  whereas 
B 

v„  (or  v„)  represents  a relaxation  rate. 
r r 
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5.  COMPUTATIONS  AND  RESULTS 

In  this  section  we  evaluate  (21)  and  (23)  for  \L  and  v-  for  given 

r r 

test  wave  parameters  (wavenumber,  frequency  and  modenumber).  The  co- 
efficients G depend  on  integrals  over  WKB  modefunctions  W . in  Appendix  A 
we  show  that  to  a good  approximation  the  modenumber  dependence  of  these 
integrals  can  be  replaced  by  a delta  function  condition  (corresponding 
to  approximate  vertical  wavenumber  conservation).  This  greatly  simplifies 
the  evaluation  of  the  decay  rates  since  use  of  the  delta  functions  reduce 
their  calculation  to  a single  sum  over  modenumber  (selects  individual 
frequency  resonance  curves)  and  a single  integral  along  each  curve. 

The  wavenumber  conservation  relations  restrict  the  allowed  region 
of  triad  interaction  to  an  open  rectangle  in  m-4  space  such  as  shown  in 
Figure  1.  The  frequency  resonance  conditions  further  restrict  allowed 
interactions  to  simple  curves  within  the  rectangle.  If  we  neglect  the 
Coriolis  frequency  f in  the  dispersion  relation  (7),  the  resonance  curves 
are  straight  lines  whose  slope  and  intercept  depend  on  modenumber  ratios. 
With  f#0  the  curves  are  bent  to  an  extent  which  depends  on  the  proximity 
to  the  m and  l axes  and  on  the  value  of  k.  Two  examples  of  resonance 
curves  are  illustrated  in  Figure  1. 

The  parameter  range  for  the  calculations  is  determined  by  the  region 

of  validity  of  the  WKB  model.  We  allow  frequencies  in  the  range  f<w<N  /3 

o 

where  the  upper  limit  may  be  varied  to  test  for  sensitivity  of  results. 

If  any  member  of  a wave  triad  has  a frequency  greater  than  the  cut-off 
value,  that  triad  does  not  contribute  to  v.  Similarly,  we  employ  a long 
wavelength  cut-off  at  100  km  so  that  the  "plane  ocean"  assumption  remains 
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valid.  For  satisfactory  convergence  (~  5%)  of  the  values  for  v it  was 
sufficient  to  retain  only  modenumbers  in  a band  of  halfwidth  25  centered 
about  the  test  wave  modenumber.  That  is,  triads  lying  on  resonance 
curves  for  which  |a-B|,  | a— y | — 25  were  included  in  the  calculations . 

Changing  the  value  of  the  frequency  cut-off  by  25%  typically  pro- 
duced less  than  a 5%  change  in  v for  frequencies  w < 15f . For  higher 

a 

frequencies , important  local  frequency  and  wavenumber  interactions 
(induced  diffusion.  Section  6)  are  blocked  by  the  cut-off,  and  the  rates 
are  diminished.  Quantitative  results  are  found  up  to  =17 f or  even  higher 
for  large  modenumbers.  Changing  the  low-wavenumber  cut-off  can  also  have 
a significant  effect  ( ~ factor  of  2),  especially  for  low  modenumbers. 

This  is  again  due  to  the  induced  diffusion  mechanism  where  one  triad 
member  has  a small  wavenumber.  The  analysis  in  Section  6 provides  a 
semi -quantitative  estimate  for  the  importance  of  this  cut-off. 

Figure  2 presents  a plot  of  the  decay  rate  Vp  against  test  wave 
frequency  for  the  GM-7 6 PSD  (10).  Each  solid  curve  is  labelled  by  the 
test  wave  modenumber  a.  The  heavy  dashed  line  represents  equality  between 
decay  rate  and  linear  frequency,  and  therefore  distinguishes  regions  of 
weak  and  strong  nonlinearity  (below  and  above  the  line,  respectively). 
Since  the  theory  assumes  a weakly  nonlinear  wavefield,  little  quantitative 
reliability  can  be  placed  on  results  above  the  dashed  line.  The  high 
frequency  "kink"  in  the  curves  is  a result  of  the  frequency  cut-off  at 
21f . 

The  results  show  some  common  trends  and  features . For  given  test 
wave  frequency,  waves  with  large  modenumber  (wavenumber)  decay  most 
rapidly:  For  given  modenumber,  the  higher  frequency  (wavenumber)  waves 


decay  fastest,  and  there  is  a low  frequency  threshold  (oj^=3f)  below  which 

the  decay  rate  decreases  very  rapidly  with  decreasing  frequency.  These 

features  are  also  exhibited  in  Figure  3 which  presents  the  sane  results  i 

as  Figure  2 except  is  plotted  as  a function  of  wavenumber.  At  high 

frequencies,  a common  (modenumber  independent)  k-dependence  for  the  decay 

rate  is  apparent. 

It  was  noted  in  Section  3 that  the  expression  (23)  for  Vp  is  not 
positive  definite:  Interactions  for  which  oj  = ±(w  -w  ) ("difference 

a 6 y 

reactions")  have  the  possibility  of  giving  rise  to  initial  growth  of  the 
test  wave  amplitude.  It  is  clear  that  if  o»  <2f,  difference  reactions  are 
the  only  type  of  interactions  possible;  therefore  vp  is  most  likely  to 
be  negative  at  these  low  frequencies.  Indeed  our  results  show  growth 
(negative  Vp)  for  small  frequencies,  especially  at  high  modenumbers. 

This  is  indicated  in  Figure  4 which  is  a contour  plot  of  the  rate  calcu- 

I 

lations.  Here  the  shaded  region  corresponds  to  growth,  the  remainder 

decay.  | 

The  calculation  of  the  fluctuation-dissipation  decay  rate  v„  is 

essentially  the  same  as  for  vp,  and  the  same  cut-off  parameters  were 

P 

used.  Numerical  estimates  of  the  correlation  time  t for  the  "noise" 

c 

term  F in  equation  (15)  show  that  t ~ w 1 so  that  the  delta  correlation 

c a 

assumption  (l8)  is  valid  in  the  same  region  as  the  weak  nonlinearity 
assumption,  \>„<u  . 

Fa  j 

In  Section  3 we  saw  that  if  \p  is  an  equilibrium  PSD  for  the  mode 

(k,a)  then  Vp(k,a)  = vp(k,a).  Similarly,  if  the  wave  field  is  nearly 

in  equilibrium,  tnen  = v„  - \j_  is  "small."  An  appropriate  measure  of 

d r r 


-# 


relative  equilibrium  is  not  vg  itself  but  rather  the  ratio 


R(k,a) 


(30) 


which  compares  the  timescale  for  evolution  of  the  spectrum  to  that  for 
relaxation  of  a single  mode. 

The  ratio  R is  plotted  against  for  selected  a in  Figure  5 using 
the  GM-76  PSD.  Apart  from  a sharp  increase  in  R at  low  frequencies  and 
at  high  modenumbers  (associated  with  vp<0),  the  small  size  of  R over 
such  a wide  range  of  and  a is  supportive  evidence  that  GM-76  is  a 
"good"  representation  of  an  equilibrium  spectrum.  Below  we  will  compare 
GM-76  with  other  spectral  forms  to  test  the  sensitivity  of  this  result. 

Olbers  [1976]  and  MoComas  and  Bretherton  [1977]  have  calculated 
energy  transfer  rates  for  various  GM  models.  They  have  also  determined 
the  direction  of  the  energy  flow  by  mapping  regions  of  positive  and  nega- 
tive Vg.  We  investigate  the  effect  of  the  variation  of  two  parameters 
in  the  PSD.  Using  GM-76  as  the  basic  form  for  ip  [see  equation  (10)], 
we  change  the  wavenumber  slope  "p"  and  the  modenumber  slope  "t . " 


A.  Wavenumber  Slope  Change 

With  the  modenumber  slope  given  by  the  reference  GM-76  value 
(t=2)  we  compute  the  ratio  R for  six  values  of  the  wavenumber  slope. 

The  results  are  presented  as  contours  in  (01^,01)  space  in  Figure  6. 
Regions  for  which  R>.lare  designated  "non-equilibrium"  and  shaded  to 
indicate  the  sign  of  Energy  flows  from  regions  of  negative  (decay) 
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to  positive  (growth)  v^.  When  R<.1,  contours  of  R are  plotted  and 
the  sign  of  vn  is  ignored.  The  area  of  this  region  and  of  the  region 

D 

Rc.Ol  indicates  how  much  of  the  spectrum  is  equilibrium. 

The  GM-76  PSD  (p=2)  is  in  equilibrium  except  for  low  modenumbers 
(a<7)  where  the  action  decays  and  for  low  frequencies  and  high  modenumbers 
(m<2.5f,  a>7)  where  the  action  grows.  Decreasing  p extends  the  v^>0 
region  to  lower  modenumbers  and  higher  frequencies;  increasing  p has 
the  opposite  effect.  Overall,  GM-76  is  closest  to  equilibrium  in  the 
high  wavenumber  region,  although  the  results  for  smaller  slopes  are  not 
very  different. 

S.  Mode-number  Slope  Change 

Here  the  procedure  is  the  same  as  above  except  that  the  wavenumber 
slope  is  fixed  at  the  C-M-76  value  (p=2),  and  the  modenumber  parameter 
t is  varied.  Figure  7 shows  the  results  in  the  same  format  as  Figure  6. 
Unlike  the  previous  case , as  t decreases  the  growth  region  is  confined 
to  higher  modenumbers  while  the  decay  region  expands.  The  GM-76  spectrum 
is  again  more  closely  in  equilibrium  for  high  wavenumbers  than  the  other 
spectra,  although  for  2<t<3  there  is  little  change  in  the  equilibrium 
region. 

Generally,  the  equilibrium  region  a>5  and  coa>3f,  is  insensitive 
to  the  slope  parameters  p and  t.  However,  these  parameters  strongly 
affect  the  division  between  regions  of  growth  and  decay  of  the  action. 

For  GM-76  our  computations  are  in  agreement  with  McComas  and  Bretherton 
indicating  that  action  flows  from  aS5  into  a >10  and  w <3f. 


-27- 


6.  SPECIAL  TRANSFER  MECHANISMS 

McComas  and  Bretherton  [1977]  have  emphasized  that  certain  classes 
of  triads  have  significant  roles  in  determining  transfer  rates.  In 
this  section  we  shall  discuss  quantitatively  the  importance  of  these 
mechanisms  for  both  the  Langevin  and  Boltzmann  rate  constants. 

We  begin  with  a discussion  of  the  mechanisms  called  induced  diffusion 
and  elastic  scattering  by  McComas  and  Bretherton.  We  discuss  these 
together  as  limiting  case  I (LCI).  As  previously,  we  let  (k,a)  be  the 
test  wave;  then  LCI  corresponds  to  the  triad  conditions 

to  (m)  = ua(k)>>Ug(J,) , 

m a k>>2.  (31) 


In  this  limiting  case  the  G coefficients  have  simple  forms,  and  we 
may  re-write  (21),  (23),  and  (2lb)  as 


vx(k,a) 


[u>2(k)-f2] 


l 

b,y 


ui2(k)-f2  A 

Y X 

<o  (k)  y(k,a) 
Y ~ 


x / d2l  <pU,$) 
K<k 


v 


2 

ka  ,2.6  ,kY 


1 <5 U (k) 
+ a 


to  (k)  ± u„(i)  ] . 
Y ° 

(32) 


The  index  "x"  is  used  to  represent  "P", 


. 


nF”,  or  "B”  where 
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N to  (k) 

Ap  = £ — 'J/(k,a), 

urOO-r 

a 


N ai  (k) 

Ap  = T,"  2 ’ 


(33) 


From  (6)  and  (9)  it  is  seen  that  Ap  is  proportional  to  wave  action. 

The  quantity  v,  , is  the  overlap  integral  of  vertical  modefunctions 

X X.K 

evaluated  in  Appendix  A. 

Equations  (32)  and  (33)  show  that  for  LCI  triads,  the  equilibrium 

condition  v_  = 0 corresponds  to 
3 


Ap  = Ap 


(3U) 


or  that  A be  independent  of  modenunber  ( equipartition  of  action). 

(For  the  case  of  induced  diffusion  this  was  noted  by  McComas  and  3vetherton 
[1977]).  The  results  (32)  and  (3*0  do  not  require  either  the  WKB  or 
Boussinesa  approximations.  The  WKB  dispersion  relation  (7)  yields  for 
(3M 

<J»(k,a)  ~ a-1 


except  near  the  inertial  frequency  which  is  precisely  the  form  of  GM-76 
(10)  in  this  domain  and  for  a>>3.  The  numerical  calculations 
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in  the  previous  section  show  equilibrium  in  this  same  region.  We  shall 
show  in  the  following  that  this  is  because  LCI  triads  dominate  the  wave 
interactions . 

The  integral  v,  appearing  in  (32)  is  evaluated  in  the  WKB 

JC-K 

approximation  in  Appendix  A.  From  this  we  see  that 


vk£k  11  ^a-B-y  + ^a+B-y  + ^a-B+y^ 


(35) 


(For  our  actual  calculations  we  use  the  more  accurate  form  given  in 
the  Appendix.)  The  terms  B = ±(a-y)  when  y=a>>B  correspond  to  induced 
diffusion  (ID).  The  term  B = ct+y  corresponds  to  elastic  scattering 
(ES). 

Analytic  evaluation  of  (32)  for  these  two  cases  gives  the  expressions 
Induced  diffusion: 
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Elastic  scattering: 


1 (a-k)ft>>a 


(36) 
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The  sun  above  represents  the  sum  over  the  two  cases 

a)  = u (k)±f  . (38) 

Y a 

For  the  analytic  evaluation  of  (36)  and  (37)  the  upper  limit  of  the 

integral  over  l in  (32)  has  been  taken  to  be  infinity  since  this  has 

little  effect  on  the  result.  The  numerical  calculations  reported  in 

Section  5 use  a lower  limit  l . = 2ir/(l00km) . This  results  in  the 

mm 

multiplication  of  expressions  (36)  and  (37)  by  a factor 

[ 1 - — tan-1( l . II  /Smf)]. 
tt  nun  o 

In  Figure  8 we  compare  the  ID  expression  (36)  (dashed  curve) 

with  the  calculated  numerical  results  (solid  curve)  for  (for  this 
comparison  the  ^-integral  cut-off  was  included).  It  is  evident  that 
except  for  f<m<3f , the  ID  mechanism  provides  a very  good  approximation 
for  Vp,  and  that  vp  may  be  easily  calculated  from  the  analytic  expression 
(36). 

In  Figure  9 we  compare  the  ES  (36)  and  ID  (39)  expressions  for 
Vp-  It  is  seen  that  ES  is  much  less  important  than  ID. 

The  parametric  subharmonic  instability  (PSI)  mechanism  of  MaComas 
and  Bvetherton  [1977]  corresponds  to  a large  scale  double  frequency  wave 


interacting  with  two  smaller  scale  waves  with  nearly  equal  frequencies. 
In  our  notation,  this  corresponds  to 


2,<<k,m 


S<<ct,y  , 


— = U)  =0) 

2 a y 


(39) 


Since  the  ID  mechanism  successfully  accounts  for  the  interactions  with 

u>  > 3f , we  expect  that,  at  most,  PSI  could  be  important  in  the  region 

a)  <2f.  For  this  case  we  find 
a 

vp  = - (1+0) 


and  an  analytic  evaluation  is  again  possible.  In  contrast  to  the  ID 
and  ES  derivations  a finite  upper  limit  must  be  imposed  on  the  £ -wavenumber 
integration  for  convergence  of  the  approximate  expressions.  With  the 
choice  l = k/5  the  PSI  limit  for  was  one  to  two  orders  of  magnitude 
less  than  the  full  computation  for  the  range  f<u  <2f.  Simulation  of  the 
PSI  conditions  (39)  in  the  full  numerical  code  has  shown  that,  although 
the  frequency  and  modenumber  limiting  conditions  were  both  well  satisfied 
by  the  important  triads , the  wavenumber  condition  was  not . We  conclude 
that  for  f<uct<2f  the  important  triads  have  comparable  wavenumbers . 


Choosing  ^max  = k,  a rather  crude  approximation  for  Vp  in  the 
range  f<wa<2f  can  be  derived  which  gives  the  correct  order  of  magnitude: 


3.5  (10)"'  k ~ ^2  [3  - tan"xQ] 

0) 

at 


ir(Ua>2-f2) 

a 


Values  calculated  from  (Ul)  are  shown  in  Figure  9. 

The  conclusions  of  this  section  are  that  for  uia>3f  anda>5,  the 
ID  mechanism,  expressed  by  (36)  provides  a good  description  of  relaxation 
within  the  internal  wave  field.  The  ES  and  PSI  mechanisms  are  never  of 
much  importance. 

Equation  (32),  its  analytic  solution,  and  the  good  agreement  with 
the  full  numerical  calculation  provide  the  explanation  for  the  conclusion 
that  the  GM-76  spectrum  is  in  equilibrium  in  the  domain  ua>3f  and 
a>5.  Outside  this  domain  the  spectrum  does  not  seem  to  be  in  equilibrium. 
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7 . ENERGY  FLOW 


Three-wave  interactions  are  only  one  of  the  many  processes  which 
contribute  to  the  overall  energy  balance  in  the  ocean.  To  understand 
the  measured  internal  wave  energy  spectrum  requires  a quantitative  measure 
of  flow  rates  between  the  various  sources  and  sinks.  Here  we  obtain  an 
estimate  of  the  energy  flow  rate  through  the  non-equilibrium  region  of 
GM-76  which  also  indicates  the  energy  requirements  to  maintain  the 
spectrum. 

The  mean  energy  per  unit  area  E(k,a)  for  the  test  wave  is  [I,  eq.  2.19] 


E(k,a)  = oj  (k)  <J,  > 
_ * a ka 


(1+2) 


The  flow  rate  is  obtained  by  taking  the  time  derivative  of  (42)  using 
equations  (24),  (6)  and  (9): 


dE(k,a) 

dt 


Z 2 ,2 
0 oj  —1 
a 


i|i(k,a)  P N^B 

D OO 


(43) 


Since  the  radiative  transfer  equation  conserves  total  energy,  the 
net  energy  flow  from  the  region  Vg<0  is  equal  to  the  flow  into  \jg>0.  We 
compute  the  flow  from  the  negative  region  by  integrating  (43)  over  this 
domain  of  k,a  space: 


dE 


(-) 


dt 


dE 

l “ 

a 


(-) 


dt 


l l 

a k 


dE(k,a) 


dt 


V° 


(44) 


The  calculated  partial  rates  dE^  Vdt  are  given  in  table  1 for 

modenumbers  between  a=l  and  a =6.  For  the  first  5 modenumbers  they  are 

comparable , while  for  a=6  the  rate  has  decreased  by  nearly  a factor  of 
2,  and  the  contributions  from  a>6  (not  shown)  are  much  smaller.  Except 
for  a=l,  the  main  contribution  to  the  k sum  in  (M)  peaks  near  u>a  - 2.2f, 
independent  of  a.  Although  the  peak  is  not  sharp  it  does  suggest  that  the 
flow  through  the  spectrum  to  smaller  frequencies  is  weak:  The  major  flow 
is  from  low  to  high  modenumbers . The  net  energy  flow  rate  from  the  region 

of  negative  v_  sums  to  dE^  Vdt  = 6.b  (10)  ^ W/m^.  for  a<6,  and  we  believe 

D 

that  this  represents  only  a slight  underestimate  of  the  total  rate. 

Others  [1976]  has  evaluated  the  transfer  rate  to  high  modenumbers 

”3  -2 

for  the  GM-75  spectrum  and  found  it  to  be  3 (10)  Wm-  , significantly 
larger  than  our  GM-76  rate.  MaComas  and  Bretherton  [1977],  however, 
have  noted  a sensitivity  of  results  to  spectral  shape  as  is  reflected 
in  our  Figures  6 and  7-  Thus  the  validity  of  computed  transfer  rates 
certainly  depends  on  the  precision  to  which  the  spectrum  is  known. 
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aFPENDIX  A 

OVERLAP  INTEGRALS  POR  EXPONENTIAL  PROFILE 

The  explicit  form  of  the  coupling  coefficients  G in  the  equations 
of  motion  (5)  is  given  in  the  appendix  of  I.  The  only  Vaisala  profile 
dependent  terms  in  these  coefficients  are  the  "overlap  integrals": 

1 ° 

v,  „ = — / pW.'  W'  W'  dz  , 

k£m  p J rr  kci  my 
o -n 


V 1 ° 

= — / pW'  w..w  dz  , 
£m  p 1 „ ka  S,3  my 
O — n 


lAl) 


where  the  modefunctions  W(z)  are  solutions  of  (1+)  with  the  normalization 


\ r 2 P 

— / (N  -r  )p  W.  W,  dz  = 6 
« J ka  k6 


o -H 


ag 


(A2) 


The  Boussinesq  approximation  entails  setting  p = Pq  = p(z=0)  in  these 
integrals  as  well  as  in  (L).  Since  typically  p varies  by  only  ~.3%  over 
the  depth  of  the  ocean,  this  approximation  has  little  effect  on  the  values 
of  the  overlap  integrals  (Al). 


If  we  restrict  ourselves  to  the  region  w<<N  the  p„  integral 

o icm 

can  he  neglected.  Qualitatively  this  can  he  seen  hy  considering  the 
constant  profile  model 


N = N 


B < z < 0 


for  which 


— - — QQQ  {6  +6„  +5  , 

N3/2B  a S ' a-8-Y  S-Y-a  Y-a-B 

o 


N3/2B 

o 


Qa  {-6a-e-Y+66-Y-a+5Y-a-B} 


From  the  form  of  the  coupling  coefficients  [I,  equation  (A.2)],  it 
seen  that  the  ratio  of  the  p to  the  v containing  terms  is  given  hy 


(k/Q)2  ~ (to/N  )2  « 1. 
o 


This  result  is  confirmed  quantitatively  by  numerical  evaluation  of  the 

overlap  integrals  for  the  exponential  profile. 

To  compute  v,_m,  it  is  sufficient  to  use  the  WKB  solutions  of  (A3) 

correct  to  first  order  in  n = u>/N  : 

o 


where 


/orT  [1-  "5“  n (k)]  -z/23 

\a(a)  = IT^a  e sin[<ji(z)  + ir/U]  , 

o a 


3sL(o) 

j .-§■ 

a 2 


Aa  = w(a-H;) , 


\aU)  = 


a N 


B , TT  N ’ 

1-  -r  n o 
2 a 


JL  JL 
N 2 n0 


i JL 
1_  2 na 


The  WKB  overlap  integral  is  formed  by  differentiating  W while  keeping 

ka 

-z/23 

e fixed,  and  substituting  into  (Al).  The  result  is 


Jkim  2^  *a  i.6  my  + ^mk  + "“du 


where 


e (-1)01  6 Y J [C  (q)  sinq  - S.(q)  cosq]  . 


Here  S2  and  C2  are  Fresnel  Integrals  [Abramouitz  and  otecun,  196k],  and 
q is  the  dimensionless  vertical  wavenumber  mismatch  at  the  surface 

q = BiV0)  -V0)  "QmY(0)l  • (A 


-Ul- 


To  first  order  in  n this  can  be  written  as 


4 = + IT  + \ (VcTVb"Vy^  ' (A10) 

The  function  I „ is  peaked  about  zero  vertical  wavenumber  mismatch. 

£m 

This  gives  an  approximate  vertical  wavenumber  conservation  law  analogous 

to  that  in  (A5)  for  the  constant  N model.  Comparing  (A7)  with  (A5) 

further  shows  that  plays  the  role  of  effective  vertical  wavenumber. 

A simple  functional  fit  to  I is  given  by 

x,m 

I^m(q)  = -i“  e"aq  cos(bq) (-l)a-e-Y  , (All) 

where  a = . 0^53  and  b = .37^3.  It  involves  an  error  of  less  than  b% 
for  0<q  <h . 

To  facilitate  the  evaluation  of  the  Langevin  decay  rates  we  approxi- 

mate  the  function  I„  with  a delta  function,  thereby  enabling  one  of  the 
im 

modenumber  summations  in  (22)  and  (25)  to  be  replaced  by  a trivial  integral. 
Using  expression  (All)  it  can  be  readily  verified  that  the  replacement 


Ik 

J'2.m 


6 (q) 


(A12) 


is  valid. 

Numerical  computation  of  the  overlap  integrals  using  the  exact 
Bessel  function  solutions  [Garrett  and  tfunk,  1972]  show  that  the  accuracy 
of  (A7)  decreases  with  increasing  frequencies  reaching  b0%  at  u)/Nq  ~ .25 
(providing  the  modenumbers  are  greater  than  one).  For  q = 0 the  u. 


-1+2- 


"_I  - .« 


2 

are  typically  less  than  20%  of  the  B v 
v decreases  rapidly  (see  All)  while  p 


; however,  as  q increases 

KXrlU 

is  relatively  constant. 


AFPENDIX  B 


ALTERNATIVE  EXPRESSIONS  FOR  RELAXATION  RATES 


Equations  (21)  and  (23)  can  be  re-written  in  a form  which  makes 
transparent  their  relationship  both  to  one  another  and  to  the  familiar 
radiative  transport  equation  of  Others  [1976],  MoComas  and  Bvetherton 
[1977],  and  Davidson  [1972]. 

To  do  this  we  use  the  hamiltonian  coupling  coefficients  r.,(k;d,m) 
defined  in  I [Appendix  A]  and  also  the  wave  action  variables  (6).  In 
terms  of  these  the  relaxation  rates  are 


<J.  > 
ka 


Vp(k,a)  = 1+ir 


6,7 


l <f 

d,m 


f (k;Jt,m) 


k-d 
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6,  . 5 (aj  +0)o-w  ) 
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x (<J„„>-<J  >)<J,  > 
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(B2 ) 


Substitution  of  these  expressions  into  v (2Ub)  yields  the  usual  expression. 
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The  transport  equation  (29)  can  also  be  written  in  terms  of  action 


variables 


dt  <Jka>  = 2vF(^’a)  [Jka  " <Jka>] 


where  Je  is  the  equilibrium  action.  Equation  (B3)  can  be  generalized 

[Others,  1976]  to  include  the  effects  of  a mesoscale  current  U(x)  by 

writing  J,  (x),  and  interpreting  -jr  as 
ka  _ ut 


d 3 • 3 . i 3 

d^=  * * si  + * ’ 3^  * 


• 3fl  ,* 

* = ST  * * = 


- — , fi  = (D  (k)  + k • u(x) 

3x  ot  . . 


FIGURE  CAPTIONS 

Figure  1 Resonance  curves  for  dispersion  relation  (7).  Dashed 
curves  are  for  f = 0. 


Figure  2 Decay  rate  versus  test  wave  frequency.  Curves  are  labelled 
by  modenumber.  Dashed  curve  represents  Vp  = w . GM-76 

is  used. 


I 

I 

•T 


Figure  3 Decay  rate  versus  wavenumber  for  GM-76. 

Figure  4 Contours  of  Vp/NQ  for  GM-76.  The  shaded  region  indicates 

negative  Vp  corresponding  to  test  wave  growth. 

Figure  5 Comparison  of  with  Vp.  The  quantity  R [equation  (30)] 
is  plotted  against  test  wave  frequency  for  selected  mode- 
number.  GM-76. 

Figure  6 Comparative  equilibria  for  different  wavenumber  slopes  p 
of  PSD  (10).  GM-76  uses  p = 2.  Non-equilibrium  regions 

(R>0.l)  are  shaded  and  the  sign  of  v„  is  given.  Otherwise 

D 

-2 

contours  of  R are  drawn:  Dashed  contours  denote  R = 10  , 

-3 

dots  denote  R = 10 

Figure  7 Comparative  equilibria  for  different  modenumber  slopes  t 
of  PSD  (10).  GM-76  uses  t = 2. 


Figure  8 


Comparison  between  the  dominant  induced  diffusion  transfer 
mechanism  [equation  (36)]  and  the  full  numerical  calculation 


of  Vp/N 


o' 


Curves  are  labelled  by  modenumber. 


Figure  9 Comparison  of  contributions  to  lvp|/N0  from  special  transfer 
mechanisms.  Analytic  expressions  for  induced  diffusion 


[equation  (36)  ],  elastic  scattering  [equation 

(37)  ] and  modified  parametric  subharmonic 


instability  [equation  (Ul)  — • — • — • ] are  evaluated. 
Curves  are  labelled  by  modenumbers. 
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